Benchmarking the variational cluster approach by means of the one-dimensional 

Bose-Hubbard model 



O 

(N 

C 

00 



•4— > 

CO 

-4-* 
03 



i 

C 

O 
o 



> 
On 

(N 

(N 
O 
O 



Michael KnapQ Enrico Arrigoni, and Wolfgang von der Linden 

Institute of Theoretical and Computational Physics, 
Graz University of Technology, 8010 Graz, Austria 
(Dated: June 21, 2010) 

Convergence properties of the variational cluster approach with respect to the variational pa- 
rameter space, cluster size, and boundary conditions of the reference system are investigated and 
discussed for bosonic many-body systems. Specifically, the variational cluster approach is applied 
to the one-dimensional Bose-Hubbard model, which exhibits a quantum phase transition from Mott 
to superfluid phase. In order to benchmark the variational cluster approach, results for the phase 
boundary delimiting the first Mott lobe are compared with essentially exact density matrix renor- 
malization group data. Furthermore, static quantities, such as the ground state energy and the 
one-particle density matrix are compared with high-order strong coupling perturbation theory re- 
sults. For reference systems with open boundary conditions the variational parameter space is 
extended by an additional variational parameter which allows for a more uniform particle density 
on the reference system and thus drastically improves the results. It turns out that the variational 
cluster approach yields accurate results with relatively low computational effort for both the phase 
boundary as well as the static properties of the one-dimensional Bose-Hubbard model, even at the 
tip of the first Mott lobe where correlation effects are most pronounced. 

PACS numbers: 64.70.Tg, 73.43.Nq, 67.85.De, 03.75.Kk 



I. INTRODUCTION 

Seminal experiments on ultracold gases of atoms 
trapped in optical lattices 1 - - — have lately turned the spot- 
light of scientific research interest on the Bose-Hubbard 
(BH) model- and its variants. The BH model exhibits a 
quantum phase transition from localized Mott phase to 
delocalizcd superfluid phase. The Mott phase is charac- 
terized by integer particle density, a gap in the single- 
particle spectral function and zero compressibility^ The 
regions in the phase diagram where the ground state of 
the BH model is in a Mott state are termed Mott lobes. 
The evaluation of the boundaries delimiting the Mott 
phase and other physical quantities is very demanding 
for the one-dimensional BH model as correlation effects 
are most important in the low-dimensional case. This is 
reflected in the special shape of the Mott lobes. Partic- 
ularly, the lobes are point shaped and a reentrance be- 
havior can be observed", in contrast to the mean field 
results 

In the present paper, we benchmark the variational 
cluster approach-^ (VCA) using the one-dimensional BH 
model. VCA is based on the self-energy functional 
approac h 14 ' 15 and is a variational extension of the clus- 
ter perturbation theor y 16 ' 17 (CPT), where the physical 
system is decomposed into clusters and the inter-cluster 
hopping is treated perturbativcly. It is crucial to inves- 
tigate the convergence properties of VCA in dependence 
of the variational parameter space, the size of the clus- 
ters and the boundary conditions used for the cluster 
Hamiltonian. The motivation for this research has been 
the increasing interest in strongly correlated bosonic sys- 
tems such as ultracold gases of atoms trapped in optical 
lattices^— and light-matter systems^ - — Lately, cluster 



methods have been benchmarked for fermionic systems 
in Ref. [2ll . where the authors used the one-dimensional 
fermionic Hubbard-Model^ which can be solved exactly 
by means of the Bethe ansatz^ in order to test the 
achievements of VCA. However, it remains an open ques- 
tion how VCA performs in the case of bosonic systems. 
Unfortunately, an exact solution of the one-dimensional 
BH model does not exist, as each lattice site can be occu- 
pied by infinitely many bosonic particles. Yet essentially 
exact density matrix renormalization group (DMRG) re- 
sults for the phase boundary delimiting the first Mott 
lobe are available^ and static properties, such as the 
ground state energy and the one-particle density matrix, 
have been evaluated using strong-coupling perturbation 
theory of high orderi 24 ' 25 In the present paper we dis- 
cuss the convergence properties of VCA for these physical 
quantities. 

The outline of this paper is as follows. In Sec. [TTJ 
the BH model is introduced. Section IIIII describes the 
most important aspects of VCA. In this section, the 
ordinary variational parameter space of the BH model 
used for cluster Hamiltonians with open boundary con- 
ditions is extended by an additional variational param- 
eter, which allows for a better distributed particle den- 
sity within the cluster and therefore drastically improves 
the results. The convergence properties of VCA for the 
phase boundary of the first Mott lobe arc investigated in 
Sec. IIVI by comparing the VCA results for distinct sets 
of variational parameters and cluster sizes with DMRG 
data from Ref. 0. Additionally, single-particle spectral 
functions and densities of states are evaluated and their 
dependence on the variational parameter space is dis- 
cussed. In Sec. |V] the ground state energy and the one- 
particle density matrix are calculated and compared with 
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strong coupling results of high order from Ref. Fi- 
nally, Sec. I VII concludes and summarizes our findings. 

II. THE BOSE-HUBBARD MODEL 

The BH Hamiltoniani is given by 

& = ~ l £ b l b i + 7 E ^ ^ - 1) - , (i) 

where b\ and &i are bosonic creation and annihilation op- 
erators at lattice site i, t is the hopping strength between 
two adjacent sites, U is the on-site repulsion and \i is the 
chemical potential, which controls the total particle num- 
ber N p = J2i 71 i = J2i K bi- The first term in Eq. (TTJ is 
considered as a sum over nearest neighbors j of site i. In 
the calculations and the forthcoming discussions we use 
the on-site repulsion U as unit of energy. 

III. THE VARIATIONAL CLUSTER 

APPROACH 

The variational cluster approach has been originally 
proposed for fermionic systems in Ref. [IH and has been 
extended to bosonic systems in Ref. [H. The main idea 
of VCA is that the grand potential il is expressed as 
a functional of the self-energy E. At the stationary 
point of the self-energy functional f2[E] Dyson's equa- 
tion for the Green's function is recovered. Unfortunately 
the functional f2[Z] cannot be evaluated directly, as it 
contains the Legendre transform of the Luttinger-Ward 
functionah 14 i 26 However, the latter just depends on the 
interaction part of the Hamiltonian and thus it can be 
eliminated by comparing S1[E] with the functional f2'[Z] 
of a simpler, exactly solvable system H' , which shares 
the interaction part with the physical system H. The 
system W is referred to as reference system. With these 
considerations one obtains for bosonic systems^ 

fi[E] = fi'[E] - Tr ^(-(G^ 1 - E)) 

+ Trln(-(G - 1 -E)), (2) 

where quantities with prime correspond to the reference 
system and Go is the noninteracting Green's function. 
The symbol Tr denotes both a summation over bosonic 
Matsubara frequencies as well as a summation over site 
indices. In order to be able to evaluate fi[E] the self- 
energy E is approximated by the self-energy of the ref- 
erence system E'. In practice this means that the func- 
tional f2[E] becomes a function of single-particle param- 
eters x of the reference system H ' 

O(x) = fi'(x) + Tr ln(-G'(x)) - Tr ln(-G(x)) . (3) 

The stationary condition on f2(x) now reads 

5§M=o. (4) 
ax 



The present formulation of VCA is not able to address 
the superfluid phase. Therefore our discussions will be 
restricted to Mott phase. A treatment of the superfluid 
phase would require an extension of the theory using 
Nambu formalism. 

In VCA the reference system H' is chosen to be a clus- 
ter decomposition of the physical system, which means 
that the total lattice of TV sites is divided into decou- 
pled clusters of size L. Formally, the decomposition can 
be achieved by introducing a superlatticc, such that the 
physical lattice is obtained when a cluster is attached 
to each site of the superlattice. The reference system de- 
fined on one such cluster is solved using the band Lanczos 
method i2L2S It can be carried out with open boundary 
conditions (obc), which is generally done in literature,— 
or with periodic boundary conditions (pbc) J£ Both cases 
are investigated in the next section. The Green's function 
of the physical system is obtained from the relation 

G~V) = G'- 1 ^) - V , (5) 

where 

V=-(/i-//)l + (T-T') (6) 

and G, G', and V are matrices in the site indices. In the 
latter relation T and T' are the hopping matrices of the 
physical and reference system, respectively. In order to 
evaluate the grand potential and the wave vector and fre- 
quency resolved Green's function G(k, ui) of the physical 
system we apply the bosonic Q-matrix formalism.— At 
first the bosonic Q-matrix formalism yields the Green's 
function of the physical system in a mixed representation, 
partly in real and partly in reciprocal space, G(k, to). 
This representation is obtained by a partial Fourier trans- 
form of Eq. (0) from superlattice site indices to wave 
vectors k of the first Brillouin zone of the superlattice. 
G(k, uj) is still a matrix in cluster site indices. Second, 
we apply the Green's function periodization proposed in 
Ref.[l6| in order to obtain the fully wave vector dependent 
Green's function G(k, u>). From the Green's function 
G(k, w) we are able to obtain the single-particle spec- 
tral function 

A(k, u) = — ImG(k, u) , (7) 
the density of states 

tfM = ^£4(k,w) (8) 

k 

and the one-particle density matrix 

C(|r 1 -r i |) = (ata i >. (9) 

The one-particle density matrix is the Fourier transform 
of the momentum distribution n(k), which can be eval- 
uated in a particularly accurate way by means of the 
Q-matrix formalism^ 
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The stationary point of f2(x) is determined numeri- 
cally by varying some or all of the single-particle param- 
eters of the reference system, see Eq. (j4}. In the case 
of the BH model the single-particle parameters are the 
hopping strength t and the chemical potential fx. Due 
to the breaking of translation symmetry introduced by 
the cluster partition, the particle density evaluated as 
a trace of G(k, tu) differs at the boundary of the cluster 
from the one inside the cluster. This is unfavorable as the 
physical system, which has pbc, should have a uniform 
particle density. However, for reference systems with obc 
this problem could be eased by adding another degree 
of freedom to the variational parameter space, which al- 
lows for a different on-site energy at the boundary of the 
cluster with respect to its bulk. Fortunately, in VCA 
the reference system H' can be extended by any single- 
particle terms, as these terms do not affect the Legcndrc 
transform of the Luttinger-Ward functional. It should be 
emphasized that due to this fact adding single-particle 
parameters to the variational parameter space does not 
affect the physical system. However, additional, physi- 
cally motivated variational parameters improve the ap- 
proximation of the self-energy Y. and therefore improved 
results for the grand potential of the physical system and 
for other physical quantities are obtained. Clearly, in the 
physical system the values of the parameters correspond- 
ing to the additionally introduced single-particle terms 
are zero, whereas in the reference system the values are 
determined by the stationary condition on the grand po- 
tential. Due to these considerations, the physical system 
is not changed when introducing a variable on-site energy 
at the boundary of the cluster. 

For the ID BH model the boundary of the cluster con- 
sists of the first and the last cluster site. After adding the 
additional boundary on-site energy the reference Hamil- 
tonian of cluster m is given by 



U' 



fi' E^ + ^'("1 + "i) 



(10) 




FIG. 1: (Color online) A system of two clusters. The inter- 
cluster hopping T is treated perturbatively. 



two cluster system is described by the Hamiltonian 



H 



H[ f 
ft H' 



The Schur decomposition of Eq. ([TU) yields 



(11) 



(12) 



All parameters of the adjacent cluster are treated on 
mean field level and thus we set H[ = (E). Next we 
have to calculate ft f , where f = b\ L b 21 + b\i b 1L . The 
notation b ma was used, where m denotes the cluster in- 
dex and a the site index within the cluster. With that 
we have 

1 1 f = (4 b 1L + b\ L b 2l) {b\ L b 21 + 4 b 1L ) 

= b 21 b XL Ah b 21 + b \ L b 21 b tl b lL 



2 bX-, 6n, b] r b-i 



bl, b 



b\ T b, 



J 21 "21 "XL "XL ~ "2X "2X "XL U 1L 

= h 2 i(2h 1L + 1) + hiL 

= n 21 (2(n 1L ) + 1) + (n 1L ) . (13) 

In the second step, we neglected the simultaneous hop- 
ping of two particles from one cluster to the other and in 
the last step we replaced the particle number operator of 
the adjacent cluster by a mean field approximation. The 
constant energy shift (tixl) can be ignored and therefore 
Eq. ([T2]) leads to 



w 

2, est 



H' 2 



2(n 



XL 



1 



(E) 



h 2 i = H' 2 + 5' h 2 i 



(14) 



where, in the second step, we replaced the fraction, which 
is just an unknown constant, by 5'. Assuming that the 
physical system has pbc and reiterating the above de- 
scribed procedure yields an extra term 5'n 2 i,. Thus we 
obtain in total 



where 5' is the additionally introduced variational pa- 
rameter. In order to retain the chemical potential at ap- 
proximately the same level as it would be without the 5' 
variation the term —j^, ™i ^ s a dded to the Hamil- 

tonian H' m as well. In VCA there is no need to justify 
the additional on-site energy in the reference system, as 
this term is not included in the physical system. How- 
ever, it is important to physically motivate the addition 
of variational parameters and, therefore, we show below 
that the additional on-site energy can be deduced from 
perturbation theory. Consider two clusters as visualized 
in Fig.[TJ where the inter-cluster hopping T from cluster 
m = 1 to cluster m — 2 is treated perturbatively. The 



H 2,est = H 2 + S'(n 21 +n 2L ) 



(15) 



In this way, we physically justified the additional on-site 
energy at the boundary of the cluster. In summary, we 
consider three possible variational parameters, namely 
the chemical potential [i, the hopping strength t and the 
additional on-site energy at the cluster boundary i5. 



IV. SPECTRAL PROPERTIES 

The first benchmark for VCA consists of a detailed 
investigation of the spectral properties of the one- 
dimensional BH model. In particular, we investigate the 
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convergence properties of VCA for the phase boundary 
delimiting the first M ott lobe with respect to distinct sets 
of variational parameters, cluster sizes, and boundary 
conditions of the reference system. Moreover, we study 
and discuss the consequences of different variational pa- 
rameters and boundary conditions of the reference sys- 
tem on the single-particle spectral functions and densities 
of states. 

In the calculations, we use the following combinations 
of variational parameters x = {fi}, {fi, t}, {[i, 6} and 
{fi,t,6}. It should be pointed out that the parameters 
of the references system are varied. Those of the phys- 
ical system are not modified. We always consider the 
chemical potential /x as a variational parameter, since it 
has been shown that the chemical potential \i must be 
varied in order to obtain the correct particle density for 
the physical system.^ . In general, CPT fails for bosonic 
systems, since the chemical potential of the reference sys- 
tems leads to erroneous densities in both, the reference 
system and the physical system accompanied by unphys- 
ical results, such as complex quasiparticle energies. 

Within the first Mott lobe, the particle density n has 
to be one. After determining the stationary point of the 
grand potential f2(x) with respect to the single-particle 
parameters x we always validate the particle density n, 
which is at T = given by 



(b) x=K t} 



N zJ E 



<2m(k) 



(16) 



k A m (k)<0 



where A m (k) are the poles of the Green's function and 
a m (k) are their spectral weights. The phase boundaries 
for the first Mott lobe are shown in Fig. [5] for refer- 
ence systems with obc, the variational parameter sets 
x = {/j,}, {/!, t}, {/!, 5}, and {fi,t, 5} and various cluster 
sizes L. The gray shaded area shown in all four subfigures 
displays DMRG results obtained from T. D. Kiihner et 
al. in Ref. 0. Additional work on the Mott to superfluid 
phase boundary can be found for instance in Refs. 0, 0, 
and [H and references therein. It can be observed that 
all sets of variational parameters yield reasonable results 
apart from the x = {/i, t} variation. The best result is 
achieved using the set x = {/i, t, 6} as variational pa- 
rameters. Particularly, the width of the phase diagram is 
approximated very well even at larger hopping strength 
t, where correlation effects are most pronounced, and the 
slope of the lobe tip is obtained correctly. At the point 
shaped lobe tip a Bcrenzinskii-Kosterlitz-Thoulcss tran- 
sition to a (quasi-long range ordered) superfluid phase 
occurs.— i^i^ The quality of the calculated phase bound- 
ary can be quantified by Xj the mean deviation of the 
VCA results from the DMRG data 



X = 



1 



E 



vX - VI 



(17) 



where pY and pf are corresponding phase boundary 
points calculated by means of VCA and DMRG, respec- 
tively, and M p is the number of phase boundary points, 
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FIG. 2: (Color online) First Mott lobe of the ID BH model 
obtained for reference systems with obc and the variational 
parameter sets (a) x = {/i}, (b) x = {/i, t}, (c) x = {/i, 5} 
and (d) x = {/j,, t, S}. The gray shaded area plotted in all 
subfigures indicates DMRG results obtained in Ref. 0. 



TABLE I: Quality y/10 -3 of the phase boundary for reference 
systems with obc and pbc shown in Figs.[2]and[7J respectively. 
The quality \ m evaluated using Eq. (|17[1 . 

L . . . number of cluster sites 
(i, t, S ... variational parameters 





obc 


pbc 


L 


/" 


n, t 


H, 5 


H, t, 5 




li, t 


2 


40.1 


32.7 


40.1 


32.7 


22.3 


32.7 


4 


23.8 


10.5 


22.7 


16.5 


16.0 


16.4 


8 


14.3 


21.1 


12.1 


8.6 


9.2 


5.8 


12 


11.5 


40.1 


8.7 


6.1 


8.0 


4.4 



which contribute to the sum. The quality \ is stated 
in Tab.U for several sets of variational parameters and 
cluster sizes. From that it can be seen as well that for 
reference systems with obc the x = {/i, t, 5} variation 
yields the best approximation for the phase boundary, as 
compared with DMRG. 

In conventional variational methods, such as Hartree 
Fock, the "best" among a set of solutions (at given /z and 
T) is chosen according to the principle of minimum grand 
potential. This criterion cannot be applied in our case, 
since there is no such minimum principle in VCA. In addi- 
tion, when including an additional variational parameter 
such as 5, there is no reason why the grand potential £1 
should display a saddle point at S — since f2 is not an 
even function of 6. 

The single-particle spectral function A(k, uj) has been 
evaluated for t = 0.15 and /i = 0.35, which corresponds 
to a point right in the middle of the first Mott lobe. It is 
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FIG. 3: (Color online) Spectral function ^4(k, u>) and density 
of states N(ui) for t = 0.15, /i = 0.35, reference systems of 
L = 12 sites with obc and the following sets of variational 
parameters (a) x = {/i}, (b) x = {/i, t}, (c) x = {/i, S}, and 
(d) x = {/i, t, 5}. 



shown in Fig. [3] along with the corresponding density of 
states N(ui) for reference systems with obc and L = 12 
sites, and distinct sets of variational parameters. For 
the numerical evaluation we used an artificial broadening 
t] = 0.05. Close to the main gap at u — /x = the spectral 
functions obtained from the {^t} and {fi, i] variation are 
not smooth but exhibit spurious gaps^ see Figs.|3](a) and 
(b). However, as soon as the variation in 6 is considered 
the spectral function becomes smooth, see Figs.|3](c) and 
(d). This also manifests in the density of states. Interest- 
ingly, in contrast to the results in one dimension, there 
are no visible spurious gaps in the spectral functions of 
the two-dimensional BH model when only the variation 
of the chemical potential x = {fi} is considered^ This 
can be explained by the fact that in two-dimensions most 
of the cluster sites are actually boundary sites as well. 
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FIG. 4: (Color online) Difference between the parameters of 
the reference system H' at the stationary point {//, t' , 5'} and 
the parameters {fi, t, 5 = 0} of the physical system H. The 
corresponding variational parameter sets are (a.*) x = {/i}, 
(b.*) x = {/i, t}, (c.*) x = {/j,, 8}, and (d.*) x = {/i, t, S}. 



The {//, t} variation has an odd behavior, see Fig.[2](b). 
Naively one would expect that the results should improve 
for larger clusters and the more variational parameters 
are used. However, this seems not to be valid for the 
variational parameter set x = t}. Therefore, this 
case has to be studied in more detail. An important as- 
pect is, that the larger the cluster the better should be 
the approximation to the physical system. Therefore, one 
would expect that the deviations of the parameters of the 
reference system H' at the stationary point {//, t', 5'} 
from the parameters {[i, t, S = 0} of the physical sys- 
tem H should decrease with increasing cluster size. At 
the same time, one would also expect the VCA results 
to come closer to the exact ones. The parameter devia- 
tions are shown in Fig.3]for various values of the hopping 
strength t and in Fig. [5] for t = 0.2 as a function of the 
cluster size L. Please notice that as a result of the prop- 
erties of the Mott phase the difference between ^ and // 
is independent of the actual value of [i provided the val- 
ues of n are restricted to the same Mott lobe. For the 
variational parameters x = {^t}, {/x, 5} and {/i, t, 5} the 
deviations decrease steadily with increasing cluster size. 
The variational parameter 5 is somewhat an exception 
as the average value of S on the cluster is SL S /L, where 
L s denotes the number of boundary sites of the cluster, 
i.e., L s — 2 for one-dimensional lattices. Correspond- 
ingly, the scaled value of S is plotted in Figs. 2] and [S] 
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FIG. 5: (Color online) Difference between the parameters of 
the reference system H at the stationary point {//, t', 5'} 
and the parameters {/j,, t — 0.2, S = 0} of the physical system 
H in dependence of the cluster size L. The corresponding 
variational parameter sets are (a) x = {/x}, (b) x = {fi, t}, (c) 
x = {fi, 5}, and (d) x = {/i, t, 5}. 



Clearly, the results for the two-site cluster are the same 
for the {/i} and {fi, 8} variation, and for the {/z, t} and 
{/^, t, 6} variation, respectively. Thus the deviations for 
that cluster size are not shown in Figs.|4](c.*) an d (d.*). 
In comparison to the deviations for all other parameter 
sets the deviations for x = {/z, t] depicted in Fig.|4](b.*) 
show a completely different behavior. In particular, they 
do not decrease with increasing cluster size and further- 
more they cross each other. This behavior suggests that 
the results might be incorrect. 

The additionally introduced variational parameter 5 
drastically improves the results and has significant influ- 
ence on the convergence properties. We suggested above 
that 5 allows for a more uniform distribution of the par- 
ticle density within the cluster. To demonstrate that this 
is indeed the case we determine both the particle density 
n'(cv) obtained from the cluster Green's function G'(oj) 
evaluated at the stationary point of the grand potential 
as well as the particle density n{a) obtained from the 
VCA Green's function G(k, uj) of the physical system. 
The latter is given partly in real and partly in recipro- 
cal space. It is important to note that at this point the 
Green's function G(k, w) has not yet been periodized, 
and that the index a is a cluster site index. Thus, it 
ranges from f to L. The particle density n{a) is evaluated 
by calculating the trace of the Green's function G(k, w), 
which reduces at T = to a sum over the residues of the 
Green's function corresponding to poles with negative en- 
ergy and a sum over the wave vectors k. The results for 
the particle densities n'(a) and n(a), respectively, for ref- 
erence systems with obc and of different size are shown 
in Fig. El As one can see from the figure, the particle 
distribution n'(a), first row, becomes flatter when <5 is 
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FIG. 6: (Color online) Particle density n'(a), first row, and 
n(a), second row, for reference systems with obc and size (a) 
L = 4, (b) L — 8, and (c) L = 12. The parameters for the 
physical system are t = 0.15 and fi = 0.35. The legend refers 
to the distinct sets of variational parameters. 



considered as variational parameter and the deviations 
from density one shrink with increasing cluster size, i. e., 
from the left to the right panel. The only exception is the 
x = {yu, £} variation, where deviations from one increase 
for larger clusters. However, a uniform particle density 
n{a) in the physical system is even more important. In 
VCA the lattice of the physical system, which has pbc 
and thus a uniform particle density, is decomposed into 
clusters of size L. This breaks the translational invari- 
ance of the physical system and hence a periodization 
prescription has to be applied "by hand" to the Green's 
function, such that it obeys the translational invariance 
of the physical lattice. The periodized Green's function 
G(k, u>) depends only on one wave vector k of the first 
Brillouin zone of the physical lattice and therefore yields, 
per construction, a uniform particle density of the phys- 
ical lattice. Nevertheless, it would be desirable to obtain 
a particle density n(a) which is as flat as possible even 
without Green's function periodization. From Fig.[6l sec- 
ond row, it can be observed that the particle density n(a) 
varies significantly with cluster lattice sites a when us- 
ing the variational parameter sets x = and {fx, t}. 
The large deviations of n(a) from 1 for these parame- 
ter sets seem to be related to the spurious gaps observed 
in the corresponding spectral functions of Figs.[3](a) and 
(b). However, as soon as 6 is introduced as variational 
parameter n(a) is indeed absolutely flat. 

Alternatively to obc we consider pbc for the reference 
Hamiltonian. The advantage of pbc is that the particle 
density n 1 (a) within the cluster is uniform. However, this 
does not necessarily imply that the particle density n(a) 
obtained from G(k, lo) is flat as well, since the additional 
hopping term between the boundary points of the cluster 
have to be subtracted again in VCA via the matrix V, 
which again breaks the translational symmetry. For ref- 
erence systems with pbc the variation in 6 appears less 
meaningful. Therefore, we consider only the variational 
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FIG. 7: (Color online) Phase boundaries of the first Mott 
lobe of the ID BH model. The reference system consists of 
clusters with pbc. The variational parameters (a) x = {/i} 
and (b) x = {[i, t} are used. The gray shaded area shows the 
DMRG results from Ref. @. 



parameter sets x = {/i} and {fx, t}. Due to the pbc there 
are additional contributions of the hopping over the clus- 
ter border which are not present in the full system. This 
contribution has to be subtracted using extra terms in 
the matrix \Z p b c , see Eq. ©. As we did throughout this 
paper, we consider here a one-dimensional lattice to de- 
duce the matrix V p f, c . However, the results can be readily 
extended to higher dimensions. The hopping matrices of 
the full system T and of the reference system T' pbc arc 
given by 



(18) 



and 



( T pbc)u - 0~pbc)manl3 — ~t 8 mn (8 a + i + 8 a +lf}) , (19) 

where {i, j} are site indices in the physical system, 
{m, n} label the clusters and {a, (3} are site indices in 
the cluster. We use pbc in the indices {i, j} and {«,/?}, 
i. e., {i, j} = N + 1 = 1 and {a, = L + 1 = 1. After 
a partial Fourier transform from superlattice site indices 
m to wave vectors k of the first Brillouin zone of the 
superlattice, one obtains 

( T - T p&cW( k ) = -(t - OOW+i + 8 a+ ip) 

+ {t-te-' l *)5 al 5p L 

+ {t-te^)8 aL 8 pi . (20) 

In contrast to obc this matrix has an extra contribution 
of — (t— t')+t = t' at its top right and bottom left corner. 
The whole matrix V p t> c (k) contains the variation of the 
chemical potential as well and hence we have 



(Vpb c ) Q/ 3(k) = -(t - t')(5 a p + i + 5 a+ ip) 
+ (t-te-^)S al S pL 



(/i - n')8 a p 



(21) 



When treating clusters with pbc by means of VCA the 
replacement of the matrix V by \Z p b c is the only new as- 
pect which has to be considered apart from the fact that 
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FIG. 8: (Color online) Spectral function A(k, w) and density 
of states N(uj) for a 12 site cluster with pbc as reference sys- 
tem. The parameters of the physical system are t = 0.15 and 
fi = 0.35. The corresponding variational parameter sets are 
(a) x = {/i} and (b) x = {fx, t}. 




FIG. 9: (Color online) Difference between the variational pa- 
rameters at the stationary point of the grand potential and 
the parameters of the physical system for reference systems 
with pbc. The variational parameter sets (a.*) x = {/j,} and 
(b.*) x = {fx, t} are used. 



the solution of the cluster itself is different due to the pbc. 
The first Mott lobe of the phase diagram for the varia- 
tional parameter sets x = {fx} and x = {/i, t}, with a ref- 
erence Hamiltonian with pbc are shown in Fig.[7J Both 
Mott lobes coincide reasonably with the DMRG data. 
The phase boundary obtained from the {fx, t} variation 
yields a very good agreement and the problems which oc- 
curred for this parameter set for reference systems with 
obc are not present anymore. The quality x of the phase 
boundary calculated according to Eq. (|T71) is stated in 
Tab.fl] From that it can be seen that the {fx, t} variation 
with pbc is comparable to the {fx, t, 5} variation with 
obc. The spectral function A(k, u>) and the correspond- 
ing density of states N(u) is shown in Fig.[S]for a 12 site 
cluster and the parameters t = 0.15 and fx = 0.35. In the 
case of pbc the spectral function and the density of states 
are not as smooth as in the case of obc and x = {fx, t, 8} 
variation but overall they exhibit the important features. 
The reason for this is that for pbc there are more conser- 
vation laws and therefore less states, so the states are less 
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FIG. 10: (Color online) Difference between the variational 
parameters at the stationary point of the grand potential and 
the parameters of the physical system for reference systems 
with pbc for the hopping strength t = 0.2 in dependence of 
the cluster size L. The variational parameter sets are (a) 
x = {fj,} and (b) x = {fj,, t}. 
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FIG. 11: (Color online) Particle density n(a) obtained from 
G(k, lo) for the parameters t = 0.15 and /i = 0.35, and refer- 
ence systems with pbc. The clusters are of size (a) L = 4, (b) 
L = 8, and (c) L = 12. The legend refers to the considered 
sets of variational parameters. 



dense. The deviation of the variational parameters with 
respect to the system parameters is as demanded shrink- 
ing for increasing cluster size, see Figs. [9] and [TO] Next 
we investigate the particle density n(a) obtained from the 
not yet periodized Green's function G(k, lu). Results are 
shown in Fig. [11] Interestingly, for a reference systems 
with pbc the density is less uniform than for a reference 
system with obc which has an additional boundary po- 
tential described by the variational parameter 5. 

Next we perform finite size scaling for the energy gap 
A L . Assuming a 1/i dependence (L is the cluster size) 
we estimate the infinite-system gap A°°. Notice that, 
in principle, VCA provides the gap for an infinite sys- 
tem. However, due to the approximation, the gap dis- 
plays a dependence on the cluster size L, which should 
converge toward the exact value for L — > oo. In Fig. [12] 
we show the deviations of the gap A°° from the gap ob- 
tained from DMRG A DMRG . One can observe that in 
this case the best results are obtained for the parameter 
set x = {/i, t, 6} and obc for the reference system. The 
scaled results for reference systems with pbc are not as 
good as the scaled results for reference systems with obc. 
At the first sight it might be counterintuitive that the 
agreement between the VCA and DMRG results is good 
close to the lobe tip. Yet, it should be noticed that the 
gap shrinks rapidly with increasing hopping strength and 
that we show in Fig. [T^] the absolute error rather than the 



FIG. 12: (Color online) Deviations between the extrapolated 
gaps A°° (VCA) and the gap A DMRG (DMRG), where (a) 
shows the deviations for reference systems with obc and (b) 
for reference systems with pbc. The dashed lines correspond 
to exact diagonalization results for the gap extrapolated to 
infinitely large systems. 



relative one. 

We also compared the extrapolated VCA results with 
extrapolated exact diagonalization (ED) results, where 
we used systems with pbc of up to 12 lattice sites. The 
best VCA results (for x = {/i, 5, i] and a reference sys- 
tem with obc) are significantly superior to the best re- 
sults obtained by extrapolating the bare ED data. The 
discrepancy is particularly pronounced deep in the Mott 
lobe for t « 0.15. 

In all spectral functions, regardless of the bound- 
ary conditions, we observe additionally to the two pro- 
nounced cosinclikc shaped bands centered around lu — fi = 
other bands with little spectral weight located at higher 
energies^ The intensity and width of these bands in- 
creases with increasing hopping strength t. In order to 
understand the additional bands we apply first-order per- 
turbation theory on the ground state, where we consider 
the hopping term of the BH Hamiltonian, see Eq. ([!} , 
as perturbation, i.e., Hi = —t^2u 7 \^l^j- The ground 



state IV'o '') °f t ne unperturbed system with particle den- 
sity n is given by 



N 



which is a tensor product of the individual single-site bare 
states | n). The energy of the bare states is Ei n \ — U n(n— 
l)/2 — /i n. The first-order correction of the ground state 
describing quantum fluctuations is of the form 



A' 



where I and I' are nearest neighbors and AE = E\ n+ i\ 



Ei 



n-l) 



2E\ n \. The first-order correction |A^)g ) is pro- 
portional to the hopping strength t, which reflects the 
fact that the intensity of the additional bands increases 
with increasing hopping strength. 
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FIG. 13: (Color online) Small black dots indicate lattice sites 
and big red dots particles, (a) single-particle term of the 
Green's function and (b) single-hole term of the Green's func- 
tion. In both situations the additional particle/hole can move 
freely on the lattice, (c) and (d), possible higher order excita- 
tions due to quantum fluctuations in the ground state of the 
single-particle term of the Green's function. 

In the first Mott lobe each lattice site is occupied 
by a single particle, provided quantum fluctuations are 
neglected. Hence, the additional particle in the parti- 
cle part of the Green's function can move freely, see 
Fig. [131(a). This yields the pronounced upper band of 
the spectral function located at uj — [i w 0.5. The 
cosinelike shape of the band is reminiscent of the dis- 
persion relation of free particles propagating on a lat- 
tice. Equivalently, the hole, introduced in the hole part of 
the Green's function, gives rise to the pronounced lower 
band, see Fig.[13](b). Considering quantum fluctuations 
in the ground state, see Eq. ([22]) . there are three energet- 
ically distinct excitations possible, namely, excitations 
from \n + 1) to |n + 2), from \n) to \n + 1) and from 
\n — 1) to |n). Next we evaluate the location of the bands 
arising from these excitations. The first investigated sit- 
uation corresponding to the excitation from \n + 1) to 
\n + 2) yields for the location of the band 

u>l = E\ n+2) + + (N — 2) E ln) - N E ln) 

n = -a« + 3E7 , (23) 

where we used n = 1 for the first Mott lobe in the second 
step. This situation, where one lattice site is occupied by 
three particles, is sketched in Fig.[13](c). The second sit- 
uation corresponding to the excitation from \n) to \n + 1) 
results in 

u 2 p = 2E [n+1) + E ln _ 1} + (N - 3) E\ n) -NE ln) 

"= 1 -/i + 2U (24) 

and is sketched in Fig.[l3](d). The third situation (ex- 
citation from \n — 1) to \n)) corresponds to Fig.fTBKa) 
and therefore contributes to the pronounced cosinclikc 
shaped band located at u — fi « 0.5. In the following, we 
compare the perturbative results for the additional bands 
with the VCA results. Particularly, we evaluate spectral 
functions and densities of states for the variational pa- 
rameters x = {(i, t 1 S}, the chemical potential fi = 0.5 
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FIG. 14: (Color online) Spectral functions, first row, and den- 
sity of states, second row, for (a) t — 0.01, fi = 0.5, (b) 
t = 0.03, (i = 0.5 and (c) t = 0.05, /i = 0.5. For the visu- 
alization we used a threshold of 0.1. The reference system 
consists of an L = 12 site cluster with obc. 



and small hopping strength t = {0.01,0.03,0.05}, see 
Fig. [HJ For small hopping strength t the perturbative 
treatment is well suited. However, the quantum fluctu- 
ations in the ground state, and thus the spectral weight 
of the additional bands as well, are small, see Eq. (|22p . 
In order to uncover the additional bands in the spec- 
tral functions we introduce a threshold for the spectral 
weight, i.e., the maximum value of the spectral function 
is restricted to this threshold and all values larger than 
or equal to this threshold are plotted in black color in 
the figures. This unveils bands with very low intensity. 
In all spectral functions shown in Fig.[T3] an additional 
band located at ojp « 2.5 can be observed. Furthermore 
a band located at uj 2 sa 1.5 is visible in Figs.fT4l(b) and 
(c). These bands match perfectly well with the pertur- 
bative results ojp and Cj 2 of Eqs. (f2"3"| and (|2~4")l . which are 
given by J)* = 2.5 and lj 2 — 1.5, respectively, where we 
employed a chemical potential (i = 0.5 and used U as 
unit of energy. 



V. STATIC PROPERTIES 

In this section, we benchmark the VCA results using 
the ground state properties of the BH model. In par- 
ticular we investigate the ground state energy Eq and 
the one-particle density matrix C(|ri — rd) = (aTo^), 
which both have been evaluated in Ref. [2J by means of 
a strong-coupling expansion of 14 th order. In VCA the 
ground state energy per lattice site is obtained from the 
grand potential f2 evaluated at the stationary point via 
the relation Eq = fi + n a*, and the one-particle density 
matrix C(|i"i — rj | ) is evaluated from the Fourier trans- 
form of the momentum distribution n(k). As for the 
spectral properties we obtain the best results for both 
the ground state energy Eq as well as the one-particle 
density matrix C(|ri — r j | ) when using the variational 
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FIG. 15: (Color online) Comparison of the ground state en- 
ergy Eq for hopping parameters corresponding to the first 
Mott lobe obtained by means of VCA (dots connected by 
lines) and strong-coupling perturbation theory (solid line). 
The lines connecting the VCA results are guidance for the 
eyes, (a) directly compares the ground state energy obtained 
from the two approaches and (b) shows their difference. 



parameters x = {/.t, t, 8} and obc for the reference sys- 
tem. Therefore, from now on we restrict the calculations 
to this set of variational parameters and obc for the refer- 
ence system, and evaluate the convergence properties of 
VCA with respect to the size of the reference system. In 
Fig.[15](a), wc compare the VCA results for the ground 
state energy Eg for hopping parameters corresponding to 
the first Mott lobe with the strong-coupling results ob- 
tained from B. Damski et al. in Ref. l24l We observe very 
good agreement between the two approaches. The devi- 
ations of the VCA results from the perturbative results 
are shown in Fig.[15](b). For the largest cluster of L = 12 
sites the deviation is even less than 0.005 for all values 
of the hopping strength t corresponding to the first Mott 
lobe. 

In Ref. the one-particle density matrix for nearest 
neighbors C(l), next nearest neighbors C(2) and next- 
next nearest neighbors C(3) have been evaluated as a 
function of the hopping strength t by means of strong- 
coupling perturbation theory and the results have been 
compared with DMRG data evaluated for an N = 40 site 
system. For C(l) the strong-coupling results compare 
well with DMRG up to the hopping strength t ~ 0.3 
and for C(2) and C(3) up to t « 0.2. In Fig.[TB] wc 
compare our VCA results with the one-particle density 
matrix of Ref. and find good agreement within the 
above mentioned regions of the hopping strength. 

The onc-particlc density matrix decays exponentially 
for large enough distances ri — r j | , i.e., 

C(|r i -r j |)oce- |ri - r j | ^, (25) 

where £ is the correlation length and its inverse describes 
the dominating exponential drop-off. The inverse of the 
correlation length l/£, shown in Fig.HBl(d). is almost zero 
when approaching the lobe tip which is already a precur- 
sor for the superfluid phase where the correlation length 
diverges. 
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FIG. 16: (Color online) One-particle density matrix C(|rj — 
rj|) for (a) nearest neighbors C(l), (b) next nearest neighbors 
C(2) and (c) next-next nearest neighbors C(3) obtained from 
VCA (dots connected by lines) and strong-coupling pertur- 
bation theory (solid line), respectively, (d) shows the dom- 
inating exponential drop-off l/£ of the one-particle density 
matrix for large enough distances in dependence of the hop- 
ping strength t and the size of the reference system. The 
inset shows an extract of l/£ where the hopping strength is 
restricted to 0.1 < t < 0.3. 



VI. CONCLUSIONS 

In the present work, we benchmarked the variational 
cluster approach by means of the one-dimensional Bosc- 
Hubbard model. In particular, we investigated the con- 
vergence properties of the variational cluster approach 
with respect to the variational parameter space, cluster 
size and boundary conditions of the reference system. 
For reference systems with obc we introduced, addition- 
ally to the chemical potential (i and the hopping strength 
t, the variational parameter 8, which allows for a modi- 
fied on-site energy at the boundary of the cluster. Using 
the additional on-site energy 8 as variational parameter 
drastically improves the results as it restores to large ex- 
tend the uniform particle density, which is violated by 
the breakup of the physical system into decoupled clus- 
ters. The resulting densities are essentially uniform, both 
that of the reference system as well as that for the physi- 
cal system computed via the variational cluster approach 
Green's function without pcriodization. 

At first we compared the variational cluster approach 
results for the phase boundary delimiting the first Mott 
lobe with essentially exact density matrix renormaliza- 
tion group data?^ Especially accurate results for the 
phase boundary have been obtained using the variational 
parameter set x = {/i, t, 8} and reference systems with 
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obc, and x = {fi, t} and reference systems with pbc. 
However, when extrapolating the results to infinitely 
large clusters the x = {/i, t, 6} variation with obc is supe- 
rior to x = {/i, t} and pbc. Naively, one would expect to 
obtain better results if additional variational parameters 
are introduced and with increasing cluster size. Both is 
not generally true. For instance, augmenting the initial 
set of parameters x = {fi} by the intra-cluster hopping 
parameter t' worsens the results. Similarly, using the 
parameter set x = {/x, t} and increasing the cluster size 
leads to monotonically increasing deviations from the ex- 
act results. This trend is accompanied by an increasing 
deviation of t' from the value t of the physical system. 
This poses the problem as to how to diagnose conver- 
gence toward the correct result, in cases where the latter 
is not known. As an indication for correct results one 
may look at the deviations of the variational parame- 
ters at the stationary point of the grand potential from 
the physical system parameters. These deviations are 
expected to shrink with increasing cluster size as larger 
clusters should better approximate the physical system. 
For x = {/i, t} and obc, however, the above considera- 
tions are not fullfilled and thus the poor results for the 
phase boundary can be ascribed to the failure of the cri- 
teria on the deviations of the variational parameters. In 
fact this criteria can be used to test the variational clus- 
ter approach results. Interestingly, for the x = {/i, t} 
variation with obc we also observe increasing deviations 
of the cluster particle density from uniform distribution 
with increasing cluster size, which might indicate that 
boundary states play an important role for this configu- 
ration. Single particle spectral functions evaluated with 
x = {/i} and {/i, t} for reference systems with obc ex- 
hibit spurious gaps.— These spurious gaps are not visible 
anymore if the additional on-site energy S is considered 
as variational parameter, which is a very important im- 
provement. Spectral functions calculated for reference 
systems with pbc are not as smooth as those evaluated 
for reference systems with obc, but overall they exhibit 
the characteristic properties. 

Second, we investigated the convergence properties of 
static quantities, such as the ground state energy and 



the one-particle density matrix. We compared our vari- 
ational cluster approach results obtained for various val- 
ues of the hopping strength t, corresponding to the first 
Mott lobe, with results obtained by means of high-order 
strong-coupling perturbation theory^ We investigated 
the convergence properties of the variational cluster ap- 
proach using the variational parameters x = {(i, t, 6} and 
reference systems with obc, since this configuration yields 
the best results for both the spectral properties as well 
as the static properties. For the ground state energy we 
found excellent agreement between the variational cluster 
approach results and the strong-coupling results. More- 
over the one-particle density matrix obtained by means of 
the variational cluster approach matches very well with 
the strong-coupling results, which are for next nearest 
neighbors and next-next nearest neighbors reliable for 
a hopping strength t < 0.2. Finally, we evaluated the 
dominating exponential decay of the one-particle density 
matrix, which is the inverse correlation length. Close to 
the lobe tip the exponential decay is almost zero, corre- 
sponding to an almost infinitely large correlation length. 
This is already a precursor for the superfluid phase. 

In summary, the variational cluster approach yields ac- 
curate results with relatively low computational effort for 
both the phase boundary and the static properties of lat- 
tice bosons in one dimension, even at the tip of the first 
Mott lobe where correlation effects are most pronounced. 
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